source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

sobj <- read_rds('/home/supervisor/mist/NaviLau/RA_zhangxuan/RA_HC_merged.rds')

sobj <- JoinLayers(sobj)

sobj <- sobj |>
  mutate(barcode = str_extract(.cell, '[A-Z]+'),
         unique_bc = str_c(barcode, '_', sample))

colnames(sobj) <- sobj$unique_bc

sobj <- sobj |>
  mutate(orig.ident = sample, barcode = NULL, unique_bc = NULL) |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

Idents(sobj) <- sobj$orig.ident

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj |> VlnPlot('nCount_RNA', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10, nCount_RNA < 40000)

sobj <- sobj |>
  quick_process_seurat()

sobj |> DotPlot(c(pbmc_markers,'TRPM2'), cols = 'RdBu', cluster.idents = T)

sobj <- sobj |>
  mutate(manual_main = case_when(
    seurat_clusters %in% c(27,14,7,17,1,6,19,18,13,23,22,16) ~ 'mono/cDC',
    seurat_clusters %in% c(21,5,15,9,8,2) ~ 'CD4T',
    seurat_clusters %in% c(24,4) ~ 'CD8T',
    seurat_clusters %in% c(3,12) ~ 'NK',
    seurat_clusters %in% c(11,10) ~ 'B cell',
    seurat_clusters %in% c(20) ~ 'Plasma cell',
    seurat_clusters %in% c(26) ~ 'Mast cell',
    seurat_clusters %in% c(25) ~ 'pDC',
    .default = 'unknown'
  ))

sobj |> as_tibble() |>
  write_source_csv('ra.pbmc.metadata')

smeta <- read_csv('mission/SLE_TRPM2_MfMo/results/ra.pbmc.metadata.csv')

sobj |> DimPlot(group.by = 'manual_main', cols = 'Paired') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  ggtitle('RA PBMC') +
  theme_jpub

publish_pdf('ra.pbmc.celltype.umap.pdf', width = 70)

sobj |> DotPlot2d('TRPM2', genotype, manual_main) +
  labs(x = 'Group', y = 'Cell type')

ra_pbmc_m2 <- read_csv('mission/SLE_TRPM2_MfMo/results/ra.pbmc.trpm2.dotplot.csv')

ra_pbmc_m2 |>
  filter(group.y != 'pDC') |>
  BubblePlot(d2 = T, size = 4) + 
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 expression') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('ra.pbmc.trpm2.dotplot')

sobj_modc <- sobj |>
  filter(manual_main == 'mono/DC')

sobj_modc <- sobj_modc |>
  quick_process_seurat(skip_norm = T)

paired20 <- c(
  "#A6CEE3", "#1F78B4",  # 浅蓝-深蓝
  "#B2DF8A", "#33A02C",  # 浅绿-深绿
  "#FB9A99", "#E31A1C",  # 浅红-深红
  "#FDBF6F", "#FF7F00",  # 浅橙-深橙
  "#CAB2D6", "#6A3D9A",  # 浅紫-深紫
  "#FFFF99", "#B15928",  # 浅黄-深棕
  "#8DD3C7", "#1B9E77",  # 浅蓝绿-深蓝绿
  "#FCCDE5", "#E7298A",  # 浅粉-深粉
  "#FFD92F", "#A6761D",  # 浅金黄-深土黄
  "#E5C494", "#8C510A"   # 浅卡其-深咖啡
)

sobj_modc |>
  DimPlot(cols = paired20, label = T, label.box = T, repel = T,
          label.size = 2) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_pdf('RA.PBMC.mono.umap.pdf', width = 60)

# Mo/MF/DC marker -----------
modc_marker <- list('cMono'=c('CD14','S100A8','S100A9'),
                    'ncMono' = c('FCGR3A','CDKN1C'),
                    'Macrophage'=c('CD68','MAFB','MARCO','CD163','MRC1'),
                    'cDC'=c('FCER1A','CLEC9A','XCR1','CD1C','CLEC10A'))

sobj_modc |>
  filter(seurat_clusters != 14) |>
  DotPlot(c(modc_marker, 'TRPM2'), cols = 'RdBu', cluster.idents = T,
          dot.scale = 2) +
  labs(x = 'Gene', y = 'Cluster') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('RA.myeloid.marker.dotplot', width = 90)

sobj_modc |>
  filter(seurat_clusters != 14) |>
  DotPlot('TRPM2', cols = 'RdBu', cluster.idents = T,
          dot.scale = 2) +
  labs(x = 'Gene', y = 'Cluster') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('RA.myeloid.TRPM2.dotplot', width = 40)

sobj_modc <- sobj_modc |>
  mutate(manual_fine = case_when(
    seurat_clusters %in% c(10) ~ 'cDC',
    seurat_clusters %in% c(19,9,5) ~ 'Non-classical Mono',
    seurat_clusters %in% c(17,2,15,6,3) ~ 'Classical Mono',
    seurat_clusters %in% c(13,1,12,4) ~ 'Classical Mo-Mac',
    .default = 'Intermediate Mono?'
  ))

sobj_modc <- sobj_modc |>
  mutate(manual_fine = case_when(
    seurat_clusters %in% c(10) ~ 'cDC',
    seurat_clusters %in% c(19,9,5) ~ 'Non-classical Mono',
    seurat_clusters %in% c(17,2,15,6,3,13,1,12,4) ~ 'Classical Mono',
    .default = 'Intermediate Mono'
  ))

sobj_modc |>
  DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('Monocyte/DC subsets') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_pdf('ra_monodc_umap.pdf', width = 70)

sobj_modc |> FeaturePlot('TRPM2', split.by = 'genotype', order = T,
                         cols = c( 'lightgrey','red'))

sobj_modc |> DotPlot(c(modc_marker, 'TRPM2'), cols = 'RdBu',
                     group.by = 'acpa_type') +
  RotatedAxis()

# TRPM2 logfc RA vs HC ------------
mono_ravhc_deg <- sobj_modc |>
  filter(manual_fine != 'cDC') |>
  FindMarkers(group.by = 'genotype', ident.1 = 'RA') |>
  as_tibble(rownames = 'gene')

mono_ravhc_deg |>
  filter(gene == 'TRPM2')

sobj_modc |>
  filter(manual_fine != 'cDC') |>
  bill.violin('TRPM2', genotype)

g1 <- last_plot()

g1 + labs(x = 'Group', y = 'Normalized expression',
          title = 'Total monocytes TRPM2') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_pdf('RA.pbmc.totalmono.M2.violin.pdf')

mono_types_ravhc <- sobj_modc |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'genotype',
                       ident.1 = 'RA')

mono_types_ravhc |>
  filter(gene == 'TRPM2')

## by sample ----------
sobj_modc |> DotPlot2d('TRPM2', orig.ident, manual_fine)

m2_expr_monotype_bysample <- last_plot() |>
  pluck('data')

m2_expr_monotype_bysample |>
  as_tibble() |>
  filter(group.y != 'cDC') |>
  mutate(group = str_extract(group.x, 'HC|RA'),
         group.y = str_remove(group.y, '\\?')) |>
  ggplot(aes(group, avg.exp, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0, size = 1) +
  facet_wrap(~group.y) +
  stat_compare_means(comparisons = list(c('HC','RA')), method = 't.test',
                     size = 2, label = 'p.signif', vjust = .05) +
  scale_y_continuous(expand = expansion(mult = c(NULL, .2))) +
  labs(x = 'Group', y = 'Average expression',
       title = 'TRPM2 in monocyte subsets') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RAvHC.pbmc.M2.expr.boxplot')

m2_expr_monotype_bysample |>
  as_tibble() |>
  filter(group.y != 'cDC') |>
  mutate(group.y = str_remove(group.y, '\\?')) |>
  left_join(acpa_meta, join_by(group.x == orig.ident)) |>
  ggplot(aes(acpa_type, avg.exp, color = acpa_type)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0, size = 1) +
  facet_wrap(~group.y) +
  #stat_compare_means(comparisons = list(c('HC','RA')), method = 't.test',
  #                  size = 2, label = 'p.signif', vjust = .05) +
  scale_y_continuous(expand = expand_scale(mult = c(NULL, .2))) +
  labs(x = 'Group', y = 'Average expression',
       title = 'TRPM2 in monocyte subsets') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('RA_ACPAvHC.pbmc.M2.expr.boxplot', width = 65)

# cell fraction change ----------
monotype_frac <- sobj_modc |>
  filter(manual_fine != 'cDC') |>
  calc_frac_conf_on_grouped_count(group = orig.ident, subtype = manual_fine)

monotype_frac |>
  mutate(group = str_extract(orig.ident, 'HC|RA'),
         manual_fine = str_remove(manual_fine, '\\?')) |>
  ggplot(aes(group, fraction*100, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0, size = 1) +
  facet_wrap(~manual_fine, scales = 'free_y') +
  stat_compare_means(comparisons = list(c('HC','RA')), method = 't.test',
                     size = 2, label = 'p.signif', vjust = .05) +
  scale_y_continuous(expand = expand_scale(mult = c(NULL, .2))) +
  labs(x = 'Group', y = '% in monocytes',
       title = 'Cell proportion') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RAvHC.pbmc.mono.subset.fraction.boxplot')

# inflam module ---------
inflam_sig_gene <- read_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sobj_modc <- sobj_modc |>
  AddModuleScore(features = list(inflam_sig_gene$gene), name = 'inflam')

# ACPA metadata ----------
acpa_meta <- read_tsv('mission/SLE_TRPM2_MfMo/data/acpa_ra_meta.tsv')

sobj_modc <- sobj_modc |>
  left_join(acpa_meta)

# save mono/DC rds ------------
sobj_modc |>
  write_rds('mission/SLE_TRPM2_MfMo/data/RA_modc.rds')

sobj_modc <-
  read_rds('mission/SLE_TRPM2_MfMo/data/RA_modc.rds')

sobj_modc |>
  DotPlot(c('inflam1', 'TRPM2'))

ra_m2_inflam <- last_plot() |>
  pluck('data') |>
  as_tibble()

ra_m2_inflam <- sobj_modc |>
  distinct(manual_fine, seurat_clusters) |>
  right_join(ra_m2_inflam, join_by(seurat_clusters == id))

ra_m2_inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = c(seurat_clusters,manual_fine)) |>
  filter(TRPM2 < .4, manual_fine != 'cDC') |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = manual_fine)) +
  geom_text_repel(aes(label = seurat_clusters), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'PBMC monocyte clusters in RA', color = 'Monocyte subtype',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RA.pbmc.mono.m2.inflam.cor', width = 70)

ra_m2_inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp,
              id_cols = c(seurat_clusters,manual_fine)) |>
  filter(str_detect(manual_fine, 'Inter'), TRPM2 < .4) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', se = F, linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = manual_fine)) +
  geom_text_repel(aes(label = seurat_clusters), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'PBMC monocyte clusters in RA', color = 'Monocyte subtype',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RA.pbmc.intMono.m2.inflam.cor', width = 70)

# M2-hi vs lo in intMono -----------
sobj_intmo <- sobj_modc |>
  filter(str_detect(manual_fine, 'Interm'))

m2hvl_intmo_deg <- sobj_intmo |>
  FindMarkers(ident.1 = 18) |>
  as_tibble(rownames = 'gene')

intmo_m2h_ravhc_deg <- sobj_intmo |>
  filter(seurat_clusters != 18) |>
  FindMarkers(group.by = 'genotype', ident.1 = 'RA') |>
  as_tibble(rownames = 'gene')

## GSEA ---------
library(clusterProfiler)

m2hvl_int_gsego <- m2hvl_intmo_deg |>
  mutate(avg_log2FC = -avg_log2FC) |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

m2hvl_int_gsego@result |>
  filter(ONTOLOGY == 'BP') |>
  select(Description) |>
  as_tibble()
  plot_enrichment(metric = NES) +
  labs(title = 'GSEA GO enrichment of RA\nTRPM2-hi intMono vs TRPM2-lo intMono') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('ra_intMono_m2hvl_gsego', width = 70)

## ORA -------------
m2hvl_int_orago <- m2hvl_intmo_deg |>
  mutate(avg_log2FC = -avg_log2FC) |>
  filter(p_val_adj < .05, avg_log2FC > 1) |>
  pull(gene) |>
  enrichGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL',
           minGSSize = 3, readable = T)

m2hvl_int_orago@result |>
  filter(ONTOLOGY == 'BP', qvalue < .05) |>
  DT::datatable()

m2hvl_int_orago <- m2hvl_int_orago |>
  simplify()

m2hvl_int_orago@result |>
  filter(ONTOLOGY == 'BP') |>
  plot_enrichment(metric = zScore) +
  labs(title = 'GO ORA enrichment of RA\nTRPM2-hi intMono vs TRPM2-lo intMono') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('ra_intMono_m2hvl_orago', width = 70)

## violin -----------
sobj_intmo |>
  mutate(cell_type = ifelse(seurat_clusters == 18,
                          'TRPM2-lo intMono', 'TRPM2-hi intMono')) |>
  bill.violin(inflam_sig_gene$gene, cell_type, facet.ncol = 5) +
  theme_pubr(base_size = 6, base_family = 'ArialMT') +
  theme(axis.text.x = element_blank()) +
  theme_jpub

publish_pdf('RA.intMono.m2hvl.inflam.gene.violin.pdf', width = 90)

# RA vs HC GSEA --------
## FindMarker across mono subset ----------
ravhc_motype_deg <- sobj_modc |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'genotype',
                       ident.1 = 'RA')

acpapvhc_motype_deg <- sobj_modc |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'acpa_type',
                       ident.1 = 'ACPA+', ident.2 = 'HC')

## cMono --------------
cMo_ravhc_gsego <- ravhc_motype_deg |>
  filter(str_detect(cluster, '^Classical'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

cMo_ravhc_gsego <- cMo_ravhc_gsego |>
  simplify()

cMo_ravhc_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  slice_sample(n = 10) |>
  plot_enrichment(metric = NES) +
  labs(title = 'GO GSEA enrichment of cMono\nRA vs HC') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('ravhc_cMono_gseago', width = 70)

cMo_ravhc_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  DT::datatable()

ravhc_motype_deg |>
  filter(gene == 'CCR2')

chemo_resp <- map_go_gene('GO:1990868')

ravhc_motype_deg |>
  filter(str_detect(cluster, '^Classic'), p_val_adj < .05, avg_log2FC > 0,
         gene %in% chemo_resp$SYMBOL) |>
  mutate(gene = fct_reorder(gene, avg_log2FC)) |>
  ggplot(aes(avg_log2FC, gene, fill = -log10(p_val_adj))) +
  geom_col() +
  labs(title = 'Top enriched chemotaxis gene:\nRA vs HC classical monocytes') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_distiller(palette = 'Reds') +
  theme_jpub

publish_source_plot('ravhc.cMono.chemotaxis.barplot', width = 70)

ravhc_cMo_ccr2 <- sobj_modc |>
  filter(str_detect(manual_fine, '^Classic')) |>
  FindMarkers(group.by = 'genotype', ident.1 = 'RA', features = 'CCR2',
              logfc.threshold = 0, min.pct = 0)

sobj_modc |>
  FindAllMarkers(group.by = 'manual_fine', features = 'CCR2')

motype_degs <- sobj_modc |>
  FindAllMarkers(group.by = 'manual_fine') |>
  as_tibble()

cMono_gsego <- motype_degs |>
  filter(str_detect(cluster, '^Classical'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

cMono_gsego <- cMono_gsego |> simplify()

cMono_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  slice_min(NES, n = 10) |>
  plot_enrichment(metric = NES) +
  labs(title = 'GO GSEA enrichment of RA cMono vs intMono+ncMono') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RA.cMono.vs.other.gsego', width = 70)

cmo_chemotax <- cMono_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0, str_detect(Description, 'chemotaxis')) |>
  pull(core_enrichment) |>
  str_split_1('\\/')

motype_degs |>
  filter(str_detect(cluster, '^Classical'), gene %in% cmo_chemotax) |>
  mutate(gene = fct_reorder(gene, avg_log2FC),
         p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  slice_max(avg_log2FC, n = 20) |>
  ggplot(aes(avg_log2FC, gene, fill = -log10(p_val_adj))) +
  geom_col() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(title = 'Top enriched chemotaxis gene:\nRA cMono vs intMono+ncMono') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('RA.cMono.vs.other.chemotaxis.barplot', width = 60)

## intMono ----------
intMo_ravhc_gsego <- ravhc_motype_deg |>
  filter(str_detect(cluster, 'Interme'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0,
        pvalueCutoff = 1)

intMo_ravhc_gsego <- intMo_ravhc_gsego |>
  simplify()

intMo_ravhc_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |> mutate(rank = rank(-NES)) |> DT::datatable()
  plot_enrichment(metric = NES) +
  labs(title = 'GO GSEA enrichment of intMono\nRA vs HC') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('ravhc_intMono_gseago', width = 70)

sobj_intmo |>
  bill.violin(inflam_sig_gene$gene, group.by = genotype, facet.ncol = 5) +
  labs(fill = 'Group', y = 'Normalized expression') +
  theme_pubr(base_size = 6, base_family = 'ArialMT') +
  theme(axis.text.x = element_blank()) +
  theme_jpub

publish_pdf('ravhc_intMono_inflam_violin.pdf', width = 90)

intMono_gsego <- motype_degs |>
  filter(str_detect(cluster, '^Interm'), p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

intMono_gsego <- intMono_gsego |> simplify()

intMono_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  slice_min(NES, n = 10) |>
  plot_enrichment(metric = NES) +
  labs(title = 'GO GSEA enrichment of RA cMono vs intMono+ncMono') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

intmo_m2h_ravhc_gsego <- intmo_m2h_ravhc_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0,
        pvalueCutoff = 1)

intmo_m2h_ravhc_gsego@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  mutate(rank = rank(-NES)) |>
  DT::datatable()

# inflam score by sample ---------
sobj_modc |> DotPlot2d('inflam1', orig.ident, manual_fine)

inflam_score_by_sample <- last_plot() |>
  pluck('data')

inflam_score_by_sample |>
  as_tibble() |>
  mutate(group = str_extract(group.x, 'HC|RA')) |>
  ggplot(aes(group, avg.exp, color = group)) +
  geom_boxplot() +
  facet_grid(~group.y) +
  stat_compare_means(comparisons = list(c('HC','RA'))) +
  labs(title = 'Inflammatory module score in RA monocyte/DC',
       y = 'Inflammatory module score') +
  theme_bw()
